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O ; Abstract 

K.^ ' We derive a stochastic nonlinear equation to describe the evolution and 

m 

QQ ' scaling properties of surfaces eroded by ion bombardment. The coefficients 

o' 

appearing in the equation can be calculated explicitly in terms of the physical 
parameters characterizing the sputtering process. We find that transitions 
C^ ' may take place between various scaling behaviors when experimental param- 

eters such as the angle of incidence of the incoming ions or their average 

penetration depth, are varied. 
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Recently there has been much interest in understanding the formation and roughening of 
nonequihbrium interfaces. Apphcations range from fluid flow in porous media to bacterial 
colonies or molecular beam epitaxy [|T| . A common feature of most rough interfaces observed 
experimentally or in discrete models is that their roughening follows simple scaling laws. 
The associated scaling exponents can be obtained using numerical simulations or stochastic 
evolution equations. The morphology and dynamics of a rough interface can be characterized 
by the surface width, U7(t, L), that scales as w^(t, L) = {[h{r,t) — h(t)]'^) = L^°/(t/L^), where 
a is the roughness exponent for the interface h{r,t) and the dynamic exponent z describes 
the scaling of the relaxation times with the system size L; h{t) is the mean height of the 
interface at time t and ( ) denotes both ensemble and space average. The scaling function 
/ has the properties f{u — > 0) ~ m^°/^ and f{u — > oo) ~ const. 

Much of the attention has focused so far on the kinetics of interfaces generated in growth 
processes. However, for a class of technologically important phenomena, such as sputter 
etching, the surface morphology evolves as a result of erosion processes 0. Motivated 
by the advances in understanding growth, recently a number of experimental studies have 
focused on the scaling properties of surfaces eroded by ion bombardment [0-|^. For graphite 
bombarded with 5 keV Ar ions, Eklund et al. ^ reported a ^ 0.2 — 0.4, and 2 ~ 1.6 — 1.8, 
values consistent with the predictions of the Kardar-Parisi-Zhang (KPZ) equation in 2+1 
dimensions P,|^. Similarly Krim et al. ^ observed a self-affine surface generated by 5 keV 
Ar bombardment of an Fe sample, with a larger exponent, a = 0.52. On the other hand, 
there exists ample evidence about the development of a periodic ripple structure in sputter 
etched surfaces (see e.g. [|]). Chason et al. have recently studied the dynamics of such 
eroded surfaces for both Si02 and Ge bombarded with Xe ions at 1 keV, and found that it 
differs from the dynamics expected for the self-affine morphologies observed in and |Q. 

In this paper we investigate the large scale properties of ion-sputtered surfaces aiming 
to understand in an unified framework the various dynamics and scaling behaviors of the 
experimentally observed surfaces. For this we derive a stochastic nonlinear equation that 
describes the time evolution of the surface height. The coefficients appearing in the equation 



are functions of the physical parameters characterizing the sputtering process. We find 
that transitions may take place between various surface morphologies as the experimental 
parameters (e.g. angle of incidence, penetration depth) are varied. Namely, at short length- 
scales the equation describes the development of a periodic ripple structure, while at larger 
length-scales the surface morphology may be either logarithmically {a = 0) or algebraically 
(a > 0) rough. Usually stochastic equations describing growth models are constructed using 
symmetry principles and conservation laws. In most cases it is difficult to derive the equation 
from a given discrete model. In contrast, here we show that for sputter eroded surfaces the 
growth equation can be derived directly from a simple model of the elementary processes 
taking place in the system. 

Most of the theoretical approaches focusing on the scaling properties of sputter rough- 
ened surfaces have assumed that essentially all relevant processes take place at the surface, 
and that nonlinear effects would appear only due to non-local effects such as shadowing . 
However, ion-sputtering is in general determined by atomic processes taking place along a 
finite penetration depth inside the bombarded material. The incoming ions penetrate the 
surface and transfer their kinetic energy to the atoms of the substrate by inducing cascades 
of collisions among the substrate atoms, or through other processes such as electronic ex- 
citations. Whereas most of the atoms which are actually sputtered are those located at 
the surface, the scattering events that might lead to sputtering take place within a certain 



layer of average depth a. As described by Sigmund's transport theory of sputtering |T^, the 
average value of the ion deposition depth depends on the energy of the bombarding ions, 
their angle of incidence, the microscopic structure of the target material and the features of 
the scattering processes taking place inside the sample. 

A convenient picture of the ion bombardment process is sketched in the inset of Fig. 
|1]. According to it the ions penetrate a distance a inside the solid before they completely 
spread out their kinetic energy with some assumed spatial distribution. An ion releasing its 
energy at point P in the solid contributes an amount of energy to the surface point O, that 
may induce the atoms in O to break their bonds and leave the surface. Following ||To| , p!T|| , we 



consider that the average energy deposited at point O due to the ion arriving at P follows 
the Gaussian distribution 

In (|I]) z' is the distance measured along the ion trajectory, and x', y' are measured in the 
plane perpendicular to it (see Fig. |l]; for simplicity in the figure x' has been set to 0); e 
denotes the total energy carried by the ion, a is the average penetration length and a and ix 
are the widths of the distribution in directions parallel and perpendicular to the incoming 
beam respectively. However, the sample is subject to an uniform flux J of bombarding ions. 
A large number of ions penetrate the solid at different points simultaneously and the velocity 
of erosion at O depends directly on the total power Eq contributed by all the ions deposited 
within the range of the distribution (||). If we ignore shadowing effects among neighboring 
points, as well as further redeposition of the eroded material, the velocity of erosion at O is 
given by 

V =p [ dr^{r) Eoir), (2) 

where the integral is taken over the region TZ of all the points at which the deposited energy 
contributes to Eq-, ^{r) is a local correction to the uniform flux J and p is a proportionality 
constant between power deposition and rate of erosion. In the following we outline the basic 



steps in the calculation of v] further details can be found in Refs. [pA] , [L2 |. 

The calculation of (0) is most efficiently perfomed in the local coordinate system (X, F, Z) 
shown in Fig. ||. The ion beam lies in the X-Z plane, forming an angle if with the Z axis, 
Z being normal to the surface at O. To simplify the calculations, we assume that: (a) 
the radii of curvature {Rxi Ry) of the surface at O are much larger than the penetration 
depth a, so that only terms up to first order in ajRxy are kept; (b) the curvatures attain 
their maximum and minimum values along either of the X and Y directions, and thus we 
can expand the value of the surface height at O by taking a second order approximation 
in the (X,y) coordinates and consistently ignoring cross term contributions. Performing 



the integral (||) the velocity v{ip, Rx, Ry) is found to be a function of the angle ip and the 
curvatures l/Rxy [0- 

Finally, from the expression of v we can obtain the equation of motion for the profile. Now 
it is convenient to use the laboratory coordinate frame (x, y, h). In the absence of overhangs 
the surface can be described by a single valued height function h{x,y,t), measured from an 
initial flat configuration which is taken to lie in the {x,y) plane. The ion beam is parallel 
to the x-h plane forming an angle < ^ < 7r/2 with the z axis. The time evolution of h is 
given by 

dh{x,y,t) / D D ^ fo\ 
— ~ -v{<^, Rx, Ry), (3) 

where ip is the angle of the beam direction with the local normal to the surface at h{x,y). 
Now (y9 is a function of the angle of incidence 6 and the values of the local slopes d^h and 
dyh, and can be expanded in powers of the latter. We will assume that the surface varies 
smoothly enough so that products of derivatives of h can be neglected for third or higher 
orders. 

At this stage additional relevant physical processes must be taken into account to describe 
the evolution of the surface. First, the bombarding ions reach the surface at random positions 
and times. We account for the stochastic arrival of ions by adding to (§) a Gaussian white 
noise 7]{x,y,t) with zero mean and variance proportional to the flux J. Second, at finite 
temperature atoms diffuse on the surface PJ^. To include this surface self-diffusion we 
allow for a term — i^V^(V^/i) [1^, where K is a temperature dependent positive coefficient. 
Expanding (|^) and adding the noise and the surface- diffusion terms we obtain the equation 
of motion |T^ 

dh dh d'^h d'^h \x ( dh\ \ ( dh\ ^^„o,„n,, 

_ . _,„ + -,_ + „ _ + u,— + - (^-j + f (^jj - A-V=(V=ft) + ,. (4) 

From (^ we can compute the expressions for the coefficients appearing in (^) in terms of 
the physical parameters characterizing the sputtering process. To simplify the discussion 



we restrict ourselves to the symmetric case a = fx. The general case is discussed in |I^. If 



we write F = (eJp/v27r) exp(— a^/2), s = sin^, c = cosO and a„ = a/a, we find for the 
coefficients in (|) 

F F 

vo = —c , 7 = —s{a\(? - 1), 
a (7 

(J ^ J 
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\y = -7-, (5) 
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Consistent with the direction of the bombarding beam and the choice of coordinates made, 
the terms in (^ are symmetric under y -^ —y but not under x —>■ —x, while for ^ ^ we get 
'J = ^x = (,y = 0, A^. = Xy and u^ = ^y |[TB|. The equation studied in Ref. |TT[] corresponds to 
the deterministic hnear version of (1), i. e. A^; = A^, = r^ = 0. 

If u^ and Uy are positive, the surface diffusion term is expected to contribute neghgibly as a 
relevant surface relaxation mechanism when we probe the system at increasingly large length 
scales. Scaling properties are then described by the anisotropic KPZ equation (AKPZ), 
which predicts two possible behaviors depending on the relative signs of the coefficients A^, 
and Xy [lT6| , p^ . If XxXy > 0, then a = 0.38 and z = 1.6, the surface width w{L,t) increases 
algebraically, being characterized by the exponents of the KPZ equation in 2+1 dimensions 
1^. For XxXy < 0, the nonlinear terms A^. and Xy become irrelevant, and the width grows 
only logarithmically, i.e. a = 0. 

In our case Uy is always negative, while z/^ can change sign as 6 and a^j are varied. A 
negative surface tension generates an instability in the system. The same kind of instability is 



known to take place in the Kuramoto-Sivashinsky (KS) equation |T8|, which is the noisyless 



and isotropic version of (1). It has been argued for the KS equation that in 1+1 dimensions 



u renormalizes to a positive value [1^, and the large length scale behavior is described by 
the KPZ equation. In 2+1 dimensions it is not completely settled whether the large distance 
behaviors of KS and KPZ fall in the same universality class, different approaches leading to 
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conflicting results pO|. 



In contrast to the KS equation, Eq. (^) is anisotropic, and explicitly contains a noise 
term. One can study Eq. (^ using a dynamic renormalization group (DRG) procedure. The 
linear instability generated by the negative surface tension leads to additional singularities 
in the Feynman diagrams of the perturbation expansion. This is related to the presence of 
a characteristic length scale in the system, generated by the competition between surface 



tension and surface diffusion, £c = v^lV'\-, where v is the largest in absolute value of 
the negative surface tension coefficients. Below we discuss the possible scaling behaviour 
predicted by (^) based on results obtained from DRG calculations performed on its isotropic 
[21[] , and the results available in the literature for other limits. The complete scaling 



version 



picture should be provided by either a DRG analysis or a numerical integration of 



The scaling behavior depends on the relative signs of i^^., z/y, A^; and A^ [^. The variations 
of these coefficients as functions of a^ and Q lead to the phase diagram shown in Fig. ^ In 
the following we discuss the various scaling regimes separately. 

Region I — For small 9 both v^ and Vy are negative. As discussed by Bradley and 
Harper [|TT|] and experimentally studied by Chason et al. IQ , a periodic structure dominates 
the surface morphology, with ripples oriented along the direction (x or y) which presents 
the largest absolute value for its surface tension coefficient. The wavelength of the ripples 
is A, = V24. 

The large length scale behavior £ ^ £c is expected to be different. Now both nonlinear- 
ities and the noise may become relevant. A DRG analysis performed for the isotropic case 
of Eqn. (^) [0,|3 suggests that inside region J, the nonlinearities X^^y flow to zero. The 
surface is logarithmically rough, with a roughness exponent a = 0. While in regions la-c 
the coefficients of the nonlinear terms may change sign as a„ or 9 vary, we do not expect 
this to modify the scaling behavior. 

Region II — This region is characterized by a positive Vx and a negative Vy. Now the 
periodic structure associated with the instability is directed along the y direction and is the 
dominant morphology at scales C, ~ Ic- Much less is known about the scaling behaviour 



for large length scales. Preliminary results of a DRG calculation indicate that the surface 
may be stabilized by the positive v^ and the nonlinearities, and (^) reduces to the AKPZ 
equation. If this is the case, then the surface roughness increases algebraically in region Ha, 
whereas logarithmic scaling characterizes region Ilh. 

Even though several aspects of the scaling behavior predicted by (^) and (^) remain to be 
clarified, we believe that these equations contain the relevant ingredients for understanding 



roughening by ion bombardment pj]. To summarize, at short length scales the morphology 
consists of a periodic structure oriented along the direction determined by the largest in 
absolute value of the negative surface tension coefficients [^. Modifying the values of a^ 
or 9 changes the orientation of the ripples [pIpT]. At large length scales we expect two 
different scaling regimes. One is characterized by KPZ exponents, which might be observed 
in region Ha in Fig. |^. Indeed, the values of the exponents reported by Eklund et al. are 
consistent within the experimental errors with the KPZ exponents in 2+1 dimensions. The 
other regime is characterized by logarithmic scaling (a = 0), which has not been observed 
experimentally so far. We believe that it can be experimentally tested for suitable (e.g, 
small) values of the angle of incidence. Moreover by tuning the values of 9 and/or a^ one 
may induce transitions among the different scaling behaviors. For example, fixing a^j and 
increasing the value of 9 would lead from logarithmic (regions la and Ic) to KPZ scaling 
{Ha), and again logarithmic scaling (Hb) for large enough angles. Also, the measurement of 
the erosion velocity may help to test the applicability of Eqn. (p. Taking spatial and noise 
averages in (^), the terms contributing to the erosion velocity are the ones proportional to 
■^0; ^x, ^y, ^x and Xy. In the stationary state, modifying 9 or a„ changes the average yield, 
allowing for a comparison of the experimental yield curves with the ones predicted by (p. 

The experimental verification of the above possibilities would constitute an important 
step to elucidate the interplay between the mechanisms leading to the different morphologies 
and dynamics for sputter-etched surfaces. It will provide as well additional insight into the 
scaling behaviors to be expected from equation (p. 
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FIGURES 
FIG. 1. Reference frames for the computation of the erosion velocity at point O. Inset: Fol- 
lowing a straigh trajectory (solid line) the ion penetrates an average distance a inside the solid 
(dotted line) after which it completely spreads out its kinetic energy. The dotted curves are equal 
energy contours. Energy released at point P contributes to erosion at O. 

FIG. 2. Phase diagram for the isotropic case o" = /U = 1. Region la: v^ < 0, A^^ < 0, A^ > 0; 
Region lb: u^ < 0, A^; > 0, \y < 0; Region Ic: Vx < 0, Xx > 0, Xy > 0; Region Ila: Vx > 0, Xx > 0, 
Aj; > 0; Region lib: Ux > 0, Xx > 0, Xy < 0. 
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